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Abstract 

We report hybrid lattice Boltzmann (HLB) simulations of the hydrodynamics of an active nematic liquid 
crystal sandwiched between confining walls with various anchoring conditions. We confirm the existence 
of a transition between a passive phase and an active phase, in which there is spontaneous flow in the steady 
state. This transition is attained for sufficiently "extensile" rods, in the case of flow-aligning liquid crystals, 
and for sufficiently "contractile" ones for flow-tumbling materials. In a quasi-lD geometry, deep in the 
active phase of flow-aligning materials, our simulations give evidence of hysteresis and history-dependent 
steady states, as well as of spontaneous banded flow. Flow-tumbling materials, in contrast, re-arrange 
themselves so that only the two boundary layers flow in steady state. Two-dimensional simulations, with 
periodic boundary conditions, show additional instabilities, with the spontaneous flow appearing as patterns 
made up of "convection rolls". These results demonstrate a remarkable richness (including dependence 
on anchoring conditions) in the steady-state phase behaviour of active materials, even in the absence of 
external forcing; they have no counterpart for passive nematics. Our HLB methodology, which combines 
lattice Boltzmann for momentum transport with a finite difference scheme for the order parameter dynamics, 
offers a robust and efficient method for probing the complex hydrodynamic behaviour of active nematics. 



1 



I. INTRODUCTION 



Active viscoelastic gels such as suspensions of active particles and active lj^u^dcry^stals^are 
soft materials receivin g in creasing theoretical and experimental attention 
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Is HVLLL8|,LL9L l2Qfl. Such materials are called "active" H21H because they 




, cell extracts H11L I12U . self-propelled colloidal parti- 
ar motors, such as actomyosin solutions 
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17LLL8D- Activity leads to striking 



continuously burn energy, for example in the form of ATP, and this drives them out of thermo- 
dynamic equilibrium even when there is no external force. Activity imparts non-trivial physical 
properties. Perhaps the most striking is that spontaneous flow can exist in non-driven active mate- 
rials , in sharp contrast to their passive liquid crystalline counterparts. Thus such 
materials, while always remaining active in a microscopic sense, can undergo a phase transition 
from a passive phase (where activity is macroscopically incoherent) to an active phase (exhibiting 
spontaneous flow). 

Active materials are typically encountered in biological contexts (although non-biological 
counterparts may also be realized, for instance with vibrated granular rods [8]). Examples include 
suspensions of bacterial swimmers 
cles 1131], and cytoskeletal gels interacting with molecu 
or microtubular networks in the presence of kinesin 11141 

phenomena such as bacterial swarming, cytoplasmic streaming and elastotaxis Furthermore, 
many biological gels, such as actin and neurofilament networks, thicken when sheared II 1911 . This 
is the opposite of the typical behaviour of viscous polymeric fluids such as molten plastics, which 
flow more easily as shear stress increases. Activity has been suggested to be amongst the possible 
causes of this peculiar flow response Jl, 20]. 

In this paper we present a series of hybrid lattice Boltzmann simulations of the hydrodynamic 
equations of motion of an active nematic liquid crystal. Derivations of the continuum equations 
we use are given in, e.g., Refs. |l], S 0] and are not repeated here. However we are aware of 
no numerical studies of the equations (with the exception of our previous work in Q20Q, which is 
a short report using a different algorithm). These are the main focus of our work. Our model 
considers a varying order parameter so that defects are automatically incorporated, as is flow- 
induced or paranematic ordering. We show that, in the limit of a uniaxial active liquid crystal with 
spatially uniform and temporally constant magnitude of order parameter (we call this limiting case 
the "Ericksen-Leslie" model in analogy with the terminology usually adopted for passive liquid 
crystals), our model reduces to the equations considered in Ref. [7]. We then consider the specific 
case of a material that is sandwiched between two infinite parallel planes at which the director field 
is anchored along a given direction. We first choose the anchoring to be along one of the directions 
in the plane (homogeneous anchoring), and we then work out the case in which there is different 
(conflicting) anchoring at the two boundary plates (homogeneous at the top, and homeotropic, i.e. 
normal to the surface, at the bottom). When the anchoring is the same at both boundaries we find 
that there is a phase transition [22] between a passive and an active phase when the "activity" 
(, a parameter which measures the coupling between pressure tensor and order parameter (see 
Section II for details), exceeds in absolute value a finite threshold. For flow-aligning materials, 
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the transition occurs for sufficiently extensile rods; for tumbling materials it occurs for sufficiently 
contractile ones. (Here "extensile" means tending to propel fluid outwards along the long axis 
or molecular director n, drawing it in radially on the midplane. while "contractile" means the 
opposite [1].) Mixed boundary conditions, instead, lead to a zero activity threshold. 

For homogeneous anchoring, we compare the numerical phase boundary to the one found in 
Ref. OD via a linear stability analysis, finding a good agreement. However, we show that the 
velocity profile found from the stability analysis is itself unstable away from the phase boundary. 

We also explore the nature of the solutions of the equations of motion (director and velocity 
field profiles) deep in the active phase, where we find that representative flow-tumbling and flow- 
aligning materials behave in a vastly different manner. The former can sustain a quasi-Poiseuille 
or banded flow, while spontaneous flow in the latter gets increasingly confined to a region close to 
the boundaries. 

Far from the phase boundary between the active and the passive phase there is strong hysteresis, 
with multistable and history-dependent solutions. These suggest that deep in the active phase the 
dynamics might be chaotic. It would be interesting to further explore the connections between 
the active nematic hydrodynamics deep in the active phase and the rheochaotic behaviour which 
selected passive liquid crystals display when they are subjected to an external forcing 11241 |25L 26]. 
There may also be qualitative analogies to the weakly turbulent viscoelastic flow discussed in [27]. 

Finally, we consider a quasi-2D case of a thin extensile flow-aligning active liquid crystal film, 
wrapped on a cylindrical surface (i.e. with periodic boundary conditions). Our simulations shows 
that there are additional instabilities in this geometry. Spontaneous flow this time appears as 
convection rolls, which, deeper in the active phase, transiently increase in number and eventually 
split up leading to a highly distorted flowing director field pattern. 

We close this introduction with some notes on nomenclature and wording. Firstly, an active gel 
is different from a fluid which is driven out of equilibrium by an external shear or heat flow, cases 
for which there is an important and vast literature (see e.g. ||28Ll29|D . In an active gel the driving 
is internal, as, for instance, a bacterium uses up ATP to propel itself. 

Secondly, at first glance our system shares some aspects with fluids which are driven out of 
equilibrium by a chemical reaction. There is a significant literature on reaction-diffusion equations 
which lead to pattern formation II23L I3CH. Bill . Ultimately, our systems are chemically driven (e.g. 
via ATP hydrolysis), but they differ from conventional reaction-diffusion systems in two ways. 
Firstly, the underlying fluid has liquid crystalline order even in the passive state. Secondly, the 
activity enters the equations of motion through a modification of the stress tensor in the Navier- 
Stokes equations by a term which is non-potential (i.e. it cannot be derived on the basis of any 
free energy). This makes the equations of active systems quite distinct from those addressed by 
reaction-diffusion models. 
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II. MODELS AND METHODS 



A. Equations of motion 

We employ a Landau-de Gennes free energy T, whose density we name /, to describe the 
equilibrium of the active liquid crystal (LC) in its passive phase (i.e. when the activity parameters 
are switched off, see below). This free energy density can be written as a sum of two terms. The 
first is a bulk contribution, 



h - ~ 7^)Qlf3 ^-QapQp-yQ-ya + "^-(Q^) 2 , (1) 



while the second is a distortion term, which we take in a (standard) one-constant approximation as 



1320 

/2 = y (d,Q a p) 2 . (2) 

In the equations above , A is a constant, 7 controls the magnitude of order (it may be viewed as an 
effective temperature or concentration for thermotropic and lyotropic liquid crystals respectively), 
while K is an elastic constant. / = f\ + / 2 is a standard free energy density to describe passive 
nematic liquid crystals I32L Here and in what follows Greek indices denote cartesian components 
and summation over repeated indices is implied. 

The anchoring of the director field on the boundary surfaces (Fig. 1) to a chosen director n° is 
ensured by adding a surface term 



h 
2 

Q% = S (n° a n°-5 a p/3) (4) 



fs = ^W (Q a p - Q%) 2 (3) 



The parameter Wq controls the strength of the anchoring, while So determines the degree of the 
surface order. If the surface order is equal to the bulk order, S should be taken equal to q, the 
order parameter in the bulk (3/2 times the largest eigenvalue of the Q tensor). W is large (strong 
anchoring) in what follows. 



The equation of motion for Q is taken to be B3I.I341. 13511 



(d t + u- V)Q-S(W,Q) = TH + AQ (5) 

where T is a collective rotational diffusion constant, and A is an activity parameter of the liquid 
crystalline gel. The form of Eq. [5] was suggested on the basis of symmetry in Refs. and 
derived starting from an underlying microscopic model in Ref. [6]. The first term on the left- 
hand side of Eq. © is the material derivative describing the usual time dependence of a quantity 
advected by a fluid with velocity u. This is generalized for rod-like molecules by a second term 

S(W, Q) = (CD + lu)(Q + 1/3) + (Q + 1/3) (fD - u) (6) 
- 2£(Q + I/3)Tr(QW) 



4 



where Tr denotes the tensorial trace, while D = (W + W T )/2 and to = (W — W T )/2 are the 
symmetric part and the anti-symmetric part respectively of the velocity gradient tensor W a p = 
dpu a . The constant £ depends on the molecular details of a given liquid crystal. The first term on 
the right-hand side of Eq. © describes the relaxation of the order parameter towards the minimum 
of the free energy. The molecular field H which provides the force for this motion is given by 

The fluid velocity, u, obeys the continuity equation and the Navier-Stokes equation, 

p(d t + updp)u a = dp(IL a p) + r)dp(d a up + dpu a ) (8) 

where p is the fluid density, 77 is an isotropic viscosity, IL a p = IT^ sslve + IT^ 1VC , and we have 
neglected an extra term proportional to d a u a which is zero in the case we are interested in (in- 
compressible fluids). The stress tensor I1^ SS1VC necessary to describe ordinary LC hydrodynamics 
is: 

passive = _ p o5af) + 2 ^ Q ^ + ±. 5oi/3 )Q ie H 1e (9) 

1 _ 1 

— iHu^Q^p + -Syp) — i{Q ai + -5 a -y)HyP 

daQyv QceyHyfi H a yQ^p 

OOpLJyv 

= CT a /3 + T a p — O a LJ yu — . 

In Eq. © we have defined the symmetric and anti-symmetric part of the passive stress tensor 
(not including the double gradient term d a Qy U Sd S g ) as a a p and r a/3 respectively, for later conve- 
nience. P is a constant in the simulations reported here. The active term is given by 

n^ tive = -CQ a p (10) 

where ( is a second activity constant . Note that with the sign convention chosen here ( > 
corresponds to extensile rods and £ < to contractile ones yj]. As for Eq. [51 the explicit form of 
the active contribution to the stress tensor entering Eq. [8]was proposed on the basis of a symmetry 
analysis of a fluid of contractile or extensile dipolar objects in yj]. It was also derived by coarse 
graining a more microscopic model for a solution of actin fibers and myosins in Ref. ||6|]. 

A full understanding of the physical origin (in both bacterial suspensions and actomyosin gels) 
of the phenomenological couplings ( and A, as well as of the range of values these may attain in 
physically relevant situations, will require multi-scale modelling at different coarse graining levels, 
and more accurate quantitative experiments. These are at the moment still lacking. However, we 
already know from experiments and from some more microscopic approaches, that actomyosin 
gels are contractile, so that in physiological conditions those materials should be described by 



negative values of ( [36]. The term proportional to A has been proposed in Ref. [1] as a symmetry 
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allowed term which, for dilute bacterial suspensions, should be negative and proportional to the 
inverse of the time scale for relaxation of activity-induced ordering. In Ref. [71] it was pointed 
out that, instead, A > when describing concentrated actomyosin gels and other systems which 
display zipping or other self-alignment effects (this is relevant for the cases considered in [37]). 

It is important to note that the model we have just written down reduces for A = £ = to 
the Beris-Edwards model for LC hydrodynamics. For a sample of uniaxial active LCs with a 
spatially uniform degree of orientational order, the director field (also called polarisation field in 
Refs. clBl H]) n is defined through 

Q a p = q(n a np - 5 a p/3) , (11) 

where q is the degree of ordering in the system (assumed to be spatially uniform). In this limit 
our model can be shown to reduce to the vectorial model considered in Jsl |5], as will be shown 
explicitly in Section III. 



B. Hybrid lattice Boltzmann algorithm 



The differential equations (0) and © may both be solved by using a lattice Boltzmann (LB) 
algorith m 113 811 . based on the 3-dimensional lattice Boltzmann algorithm for conventional liquid 
crystals |4oll. generalised to include the two extra active terms, as we discussed in Ref. |2o|l . 

Here we use a different route, and solve Eq. © via a finite difference predictor-corrector 
algorithm, while lattice Boltzmann is used to solve the Navier-Stokes equation, ©. With respect 
to a full LB approach [3^, 40], the primary advantage of this method is that it will allow simulations 
of larger systems as it involves consistently smaller memory requirements. Indeed, while in a full 
LB treatment one has to store 6 sets of 15 distribution functions at any lattice point (if we choose 
the 3DQ15 velocity vector lattice 113 811 as we do here), just one set of distribution functions plus 
the five independent components of the Q tensor, is needed in this hybrid algorithm. Furthermore, 
we avoid in this way the error term arising in the Chapman-Enskog expansion used to connect the 
LB model to the order parameter evolution equation in the continuum limit 13911 . 

Lattice Boltzmann algorithms to solve the Navier-Stokes equations of a simple fluid are defined 
in terms of a single set of partial distribution functions, the scalars fAx), that sum on each lattice 
site x to give the density. Each /$ is associated with a lattice vector e*i |40j]. We choose a 15-velocity 
model on the cubic lattice with lattice vectors: 



-(2) 



(0,0,0) 

(±i,o,o),(o,±i,o),(o,o,±i; 
(±i,±i,±i). 

corresponds to ef% i = 1 



(12) 
(13) 
(14) 

The indices, i, are ordered so that i = corresponds to ef\ i = 1, • • • , 6 correspond to the ef 3 

_/2) 

set and i = 7, • • • , 14 to the e\ set. For our hybrid code, the input to the equilibrium distribution 
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functions has to come from the solution (via finite difference methods) of the coupled Eq. Q. 
This differs from the fully LB treatment of nematics; see Refs. 139.1401. 
Physical variables are defined as moments of the distribution functions: 

i i 

The distribution functions evolve in a time step At according to 

fi{x + e t At, t + At)- fi(x, t) = Y [C f i(x, t, {fi}) + C fi (x + %At, t + At, {/*})] . (16) 

This represents free streaming with velocity e*j followed by a collision step which allows 
the distributions to relax towards equilibrium. The /*'s are first order approximations to 
fi(x + 6iAt,t + At), and they are obtained by using AtCfi(x,t, {fi}) on the right hand side 
of Eq. (fT6l) . Discretizing in this way, which is similar to a predictor-corrector scheme, has the 
advantages that lattice viscosity terms are eliminated to second order and that the stability of the 
scheme is improved [39]. 

The collision operators are taken to have the form of a single relaxation time Boltzmann equa- 
tion, together with a forcing term 

C fi (x,t,{fi}) = --{f i (x,t)- f^(x,t, {fi})) + Pi (x,t, {fi}), (17) 
T f 

The form of the equations of motion follow from the choice of the moments of the equilibrium 
distributions ff 9 and the driving terms Moreover, fl q is constrained by 

£ fi 9 = P, 2 /r e i« = P U », fi qe ia e iP = ~^aP + pU a Up (18) 



where the zeroth and first moments are chosen to impose conservation of mass and momen- 
tum. The second moment of f eq is determined by a a p, whereas the divergences of r a/ 3 and of 



d a Qju xJn enter effectively as a body force: 



Vpi = 0, Vpje ia = dpTap - dp (d a Q JU ^ ) , Vpie ia e i/3 = 0. (19) 

Conditions (fT8l) -([T9l) are satisfied by writing the equilibrium distribution functions and forcing 
terms as polynomial expansions in the velocity. The coefficients in the expansion are (in general 
non-uniquely) determined by the requirements that these constraints are fulfilled (see Ref. 1I40I1 
for details). The active contributions then simply alter the constraints on the second moment of 
the /j's. (Alternatively the derivative of the active term could be entered as a body force and thus 
would modify the constraint on the first moment of the p/s; we do not pursue this here.) 

In Appendix (0 we give a quantitative comparison between the hybrid LB algorithm used 
here and two versions of a fully LB-based code for active nematics [|20j]. The hybrid code is 
quite satisfactory in performance; it is also easier to code and runs substantially faster due to the 
elimination of the cumbersome additional distribution functions required to represent the order 
parameter dynamics within a fully LB-based approach. 
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III. MAPPING TO ERICKSEN-LESLIE LEVEL EQUATIONS 



In this section, we consider the limit of the equations of motion © and © when the ac- 
tive molecules are uniaxial, so that the order parameter can be written in the form Q af3 = 
q (n a np — 5 a p/3) (n being the usual nematic director field). We furthermore assume that the 
magnitude q of the nematic ordering is independent of space and time. The resulting simplified 
theory is commonly employed in the physics of active gels (see e.g. Ref. Jsllsl]); using it, some 
analytical results have been found. It is thus useful to explicitly consider this limit (i) to show 
that our equations map onto those of Ref. |0, Q] for uniaxial systems, and (ii) to quantitatively 
check our numerical results against those found analytically for the phase boundaries separating 
the active and passive states tn\- hi this Section quantities labelled by "EL" refer to the resulting 



director- field model, which is the direct counterpart of the Ericksen-Leslie theory B2I1 of passive 
liquid crystal hydrodynamics. 

A. Order parameter equation of motion 

We first note that the evolution equation © of the tensor order parameter can be written in the 
usual form for a purely passive system 

(d t + u-V)Q- S(W, Q) = rtT (20) 

so long as we write an effective molecular field 

H' = H + ^Q (21) 
This implies that the the classical linear (in Q) term of the molecular field, namely 

-A (l- 7 /3)Q Q/3 , (22) 



is now effectively replaced by 



f -A (1 - 7 /3) + ±\ Q a p. (23) 



In this manner the "equilibrium" properties of active nematics can be said to differ from the passive 
ones because of the presence of the active parameter A. (This contrasts with the role of £, which 
has no equilibrium counterpart. We will see below, moreover, that the shift created by A has no 
dynamical consequences in systems where the ordering strength q is fixed.) 

After some straightforward algebra (see e.g. Q33Q and Refs. therein), one finds that the linear 
term now changes sign for 7 = 7*, with 



7* = 3 1 - — - . (24) 
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Similarly the transition point 7 = 7 C for the first-order isotropic-to-nematic transition obeys 

^^-r^HK'-fy- <25) 

Furthermore, for uniaxial nematics with a spatially uniform degree of ordering q (as assumed at 
Ericksen-Leslie level - see above), the solution for q becomes 




i( x ) = i + i\ ^-tt + tt^t- ( 26 ) 



(Note that this is 3/2 times the largest eigenvalue of the Q tensor.) 

The conventional passive case is recovered by setting A = in equations (|23T[26l) . Note that 
the value of q at the transition, q c = 1/3, is independent of A since it is insensitive to the quadratic 
term of the free energy density. However, the condition for real solutions (positivity of the term 
inside the square root) for active nematics is shifted by nonzero A and becomes 

The dynamics of the director field in a uniaxial active liquid crystal of fixed q is controlled 
by three parameters. These are jel, the liquid crystal rotational viscosity; uel, which is an- 
other viscosity determining whether the liquid crystal (in its passive phase) is flow-aligning or 
flow-tumbling (for \vel\ larger and smaller than 1 respectively); and \el, which determines the 
magnitude of activity-induced ordering. It is possible to map the dynamical equation of motion 
for Q a/3 (l20l) onto the model considered in Ref. [7] (the details are worked out in Appendix |B]), 
which leads to the following identifications: 

2q 2 

Jel = 7l = -p-, (28) 

„„ = 2L = Jl±M, (29) 
7i 3g 

Xel = 0. (30) 

These relations show that in our model the dynamics of the tensorial order parameter may be 
controlled by tuning £ and T. Furthermore, we note that our parameter A does not control \el 
directly, because in Ref. [7] this parameter can already be adsorbed into a Lagrange multiplier 
introduced to maintain fixed q. (To emphasize this, we set it to zero above; see also Appendix|B]) 
However, changing A in our equations does alter q, so qualitatively the meaning of this parameter is 
similar to that of Xel m 0] insofar as it determines the strength of activity-induced self-alignment 
effects. The relations (|28l) , (|29l ), give rise to a non-trivial dependence of the parameter i/ EL on 7 
and £ as shown in figure [TJ 
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B. Navier-Stokes equation 



We now map out the parameters entering the Navier-Stokes equation ([8]) onto the analogous 
equation derived at director-field level in Ref. [0], which is written in in terms of the "vectorial" 
molecular field and of the director field, n^. In Ref. [3], the velocity field at steady state of an 
active gel is determined by pel (see SeclQ, t)el, which is an isotropic viscosity similar to the one 
introduced in Eq. ([8]), and (el, which controls the hydrodynamics in the active phase, determining 
whether the active liquid crystal is extensile or contractile as discussed in Section II. (Note that 
(el controls the effect of activity on the Navier-Stokes sector, but does not enter directly the order 
parameter dynamics as set up in SecjO) 

After some algebra (the details of which are worked out in Appendix B), we can rewrite Eq.® 
in the required limit of uniaxiality and fixed q. We find that the six Leslie viscosities for a purely 



passive liquid crystal (A = ( = 0), which are usually called cui,...^ IB2I1 . are: 

oci = -^g 2 (3 + 4g-4g 2 )e 2 , (31) 

012 = f(-^( 2 + ( ?)£-? 2 )> < 32) 

as = ^(-^(2 + g)e + g 2 ), (33) 

"4 = ^(l-q) 2 e + V, (34) 

« 5 = ^0?(4-g)£ 2 + g(2 + g)O, (35) 

«e = ^(<?(4-g)£ 2 -g(2 + g)0. (36) 

The Parodi relations, 

«3 - «2 = ~Y = 71 > ( 37 ) 

a 6 - a 5 = _ ) = 72 , (38) 

a 2 + a 3 = a 6 - a 5 , (39) 

are easily seen to hold. The Ericksen-Leslie level viscosity and active stress term are recovered as: 

vel = v + ^(q-i) 2 e, (40) 

(el = (■ (41) 



Using the above relations and the results of Ref. 07Q, we obtain the phase boundary in the ((, A) 
plane, for an active nematic confined between parallel plates at separation L, with homogenous 
anchoring at the walls (Figj3]): 

^ r2 _ U^KjUrfT -5^ 2 -14£g + e + eV + 4e 2 + 4e 2 g + 9g 2 ) 

C 9(£g + 2£-3g) " ^ 
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From Eq. (l42l) it is apparent that the critical activity threshold beyond which spontaneous flow is 
found scales like L~ 2 , and thus vanishes for an infinite system. Note that the dependence on A of 
the phase boundary is indirect, via q. Fig. [2] shows an example of comparison between analytical 
and simulated phase boundary, from which it is apparent that there is a good agreement. 

IV. RESULTS 

Most of the results which we present below refer to a quasi- ID system in which the active 
nematic is sandwiched between two plates at separation L in the z direction, with translational 
invariance assumed in x and y (Fig. [3]). We consider two different boundary conditions: either 
homogeneous anchoring along the y— direction, or mixed (conflicting) anchoring at the two plates. 
We will also refer to the angle between the director field and the positive y direction as the polar- 
ization angle, 6, the convention being that 9 > if the positive y axis can be superimposed with 
the director field with an anti-clockwise rotation of an angle 16*1 (which is defined to be smaller 
than 7r), around the x axis. 

A. Spontaneous flow transition in Freedericksz cells 

We first consider homogeneous anchoring where the polarization at the confining surface is 
parallel to the y-direction, 9 = 0. (This geometry is known as the Freedericksz cell in passive 
liquid crystal device terminology, 03211 .) By considering Eq. [26] we see that the order parameter 
q remains between and 1 for small values of A. Furthermore, we note that for £ = 0.7 and 
£ = 0.5 the system is respectively in the flow-aligning regime (point A in Fig. [U and in the 
flow-tumbling regime (point B). Let us first concentrate on the flow-aligning regime (point A). 
For definiteness we now fix A = 0, r/ = 2.5, A = 0.1, K = 0.04, V ~ 0.34 and 7 = 3; while 
£ can take on the discrete values 0.5, 0.7 as just described, and L and £ are variable. Note that, 
as described previously, setting A = eliminates the shift in q arising from self-alignment but 
this term can anyway be adsorbed into an effective (quasi-passive) free energy. Accordingly, the 
important activity parameter, for our purposes, is simply £. 

1. Flow-aligning regime 

For £ = 0.7, the system is flow-aligning and, for instance with ( = 0.005, the active LC is 
extensile. In Figure |4] we show the time evolution of the components n y , n z of the polarization 
vector at the center of a system of size L = 100 lattice units (n x is identically zero in this case). 
The polarisation field was inizialized along the y direction except for the midpoint director field, 
which was initialised with 9 = 10°. As one can see for t > t* ~ 10 5 timesteps, the system 
undergoes a transition to an active state, characterised by a spontaneous flow. 
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This happens when the scaling variable (L 2 becomes larger than the critical value found 
through the solution of Eq. |42l Thus there are two ways of entering the active phase: either 
by increasing the value of £ at fixed L, or by increasing the system size at fixed activity. In Figs. [5] 
and [6] we explore the system behaviour (respectively director and flow field at steady state) when 
the active phase is entered via an increase in the activity parameter £. 

By means of a stability analysis, valid very close to the phase boundary, an analytic expression 
for u y (z) was found in [7]. This predicts a sinusoidal modulation with a node at the centre of 
the channel. While our numerics shows this solution to be metastable for a long time close to 
the threshold, the eventual steady state we find is a quasi-Poiseuille flow with a maximum flow 
velocity, not a nodal point, at the centre of the channel (Fig. |5]). Thus with homogeneous boundary 
conditions and assumed translational invariance along the flow direction, we obtain a spontaneous 
net mass flux rather than the balancing fluxes of forward and backward fluid in the two halves of the 
cell, suggested by the analysis of [7]. Our numerical simulations thus suggest that the perturbative 
solution is stable at most within a very narrow region close to the phase boundary. The overall mass 
flux is set in a direction chosen by spontaneous symmetry breaking or, in practice, small deviations 
from symmetry between y and — y in the initial condition. Note that for a fixed initial condition 
as selected above, the flow direction can also switch on variation in (: to ease comparisons, some 
such switches are silently reversed in the figures presented here and below. 

Upon increasing the value of (L 2 (i.e. moving deeper inside the active phase) the flow pattern 
changes from quasi-Poiseuille flow to a "banded" flow, with regions of rather well defined and 
distinct local shear rates (Fig. [6]). These bands (which are clearer and more numerous in larger 
samples, see Fig. |7J correspond to regions of aligned liquid crystal, which are separated by sharp 
interfaces. As the equations deep in the active phase are strongly non-linear, no analytical results 
so far exist to probe the behaviour of an active gel in this regime. The utility of a robust numerical 
algorithm, as we have developed here with our HLB code, is highly apparent when addressing the 
potentially complex behaviour in such regimes. The model we consider allows for a non-constant 
value of the order parameter q and we can thus quantify the variations in q that are neglected in a 
director field model. Variations in q are at most of 1 — 5% in the simulations reported above, and 
small dips in the order parameter correspond to the spatially rapidly varying regions in the director 
field profile (i.e. in the "kinks" which appear at the band edges). Furthermore, these small changes 
are only encountered far from the phase boundary. 

2. Flow-tumbling regime 

We now turn our attention to the flow tumbling regime by considering £ = 0.5, 7 = 3 and 
A = 0. In this case Eq. (|42l) suggests that, in order to have a spontaneous flow, £ must be 
negative (i.e. the LC has to be contractile). This is confirmed by our simulations. We consider 
the value £ = —0.0025, which is just in the active phase (see Eq. (l42~l)). In Figure [8] we show 
the time evolution of the components n y , n z of the polarization vector at the center of a system of 
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size L = 100, inizialized as for the flow aligning case. As in the flow-aligning case, for t > t* 
the system undergoes a spontaneous alignment with a consequent spontaneous flow. The time 
behavior is however quite different from the one observed in the flow aligning case. In particular 
at t — t* the polarization vector has an abrupt variation of n/2 and then reaches a stationary value 
with a polarization angle which strongly deviates from the starting configuration. 

As with the flow-aligning case, we can estimate the critical value ( c at fixed L (or L c at a given 
above which the system starts to display spontaneous flow in steady state. Again as in the flow- 
aligning case we find good agreement between the value of the threshold estimated numerically 
and the analytical prediction of Eq. (l42l) . However, a comparison between the stationary profile 
of velocity and polarization angle profile in the flow-aligning regime and in the flow-tumbling one 
(Figs. [5] and [9] respectively) shows a striking difference. While the velocity profile has the shape 
of a spontaneous Poiseuille flow for a flow-aligning active liquid crystal, it is zero in the centre 
of the channel and confined to the boundaries in the flow-tumbling case. Also the polarization 
angle is quite different: in the flow-aligning case the director field splays and bends so that the 
polarization angle approaches the Leslie values (selected by the local shear), while it is almost 
constant throughout the sample in the flow-tumbling case. 

Upon moving deeper inside the active phase, first the velocity field becomes confined more and 
more to the boundaries, while the polarisation angle becomes increasingly close to 90° throughout 
(Fig. HI). For still larger values of the activity parameter ( (Fig. [TOl) . the flow changes sign, passing 
through an intermediate state with plug-like flow in which the polarization has the shape of a kink 
(notice however that 9 = ±90° are equivalent due to the head-tail symmetry of the director field). 
As in the flow-aligning case, order parameter variations are limited for ( just larger (in absolute 
value) than the critical value. For the simulations presented here and deep in the active phase, 
the order parameter shows some drops (similar in magnitude to those found with flow-aligning 
materials) close to the boundary plates, where the shear rates are maximal. 

3. Multi-stability in the active phase 

It is important to consider whether the solutions we have found are unique (modulo the trivial 
bistability associated with sign-reversal, discussed above), or whether each of them is one of many 
possible solutions of the equations of motion with given anchoring conditions at the boundary. The 
selection between such solutions, if they exist, is presumably governed by the initial conditions. 
We focus here, for definiteness, on the case of contractile active tumbling liquid crystals. 

Figs. [TT] and [12] show the results of two different initial conditions on the steady state director 
and velocity profiles. Fig. [TT] shows data for a modest value of the activity (~ 50% larger in 
absolute value than the critical value to enter the active phase). It can be seen that one of the 
solutions has a non-zero component of the director field along the x direction, so that the director 
tilts out of the "shear plane" (the yz plane in Fig. [3]). Fig. [12] shows another example, deeper in 
the active phase, in which the polarisation profiles again differ in steady state for the two different 
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initial conditions. One of these initial conditions is the same as above, for the other we started 
the director field along the z direction apart from (the boundary and) the midplane in which the 
polarisation angle was tilted. 

Extensile aligning liquid crystals behave in a similar way. As a rule of thumb, multi stability 
appears to increase for intermediate values of the activity. For the cases considered here, we only 
find a single (bistable) solution in the active phase close to the phase boundary and again for very 
large activity. It should be noted that also passive liquid crystals can have metastable multiple 
solution in equilibrium (for instance super-twisted structure are metastable). However, in that 
case (in the presence of thermal noise, and in the absence of external driving) one can speak of a 
"most stable solution" which is unambiguously determined by free energy minimization. No such 
criterion exists for our non-equilibrium problem, as the equations of motion cannot be written 
down completely in terms of a free energy. (Note however that, were C, = 0, this could be done 
even in the presence of the active self-alignment term A.) 



B. Spontaneous flow in hybrid aligned nematic cells 

^ow we consider a hybrid-aligned nematic cell (HAN cell, in passive liquid crystal terminology 
114 ill ), in which the polarization vector is anchored homogeneously at z = and homeotropically 
at z = L. We restrict attention to £ = 0.7, the flow-aligning case. 

Unlike the Freedericksz cell, the conflicting anchoring now leads to an elastic distortion in 
equilibrium even within the passive phase of the active system (as it would in a strictly passive ne- 
matic). As a result any non-zero value of (, whether positive or negative, leads to spontaneous flow 
in steady state, as the active pressure tensor is no longer divergence-free when C 7^ 0. Thus even 
contractile aligning liquid crystals flow spontaneously in this geometry (Fig. [13]). The velocity 
profiles in steady state in this case show extended regions with very low shear rate and plug-like 
flow, coexisting with strongly sheared "boundary layers". This is similar to what was observed 
in Section IV A.l for contractile (tumbling) liquid crystals in a Freedericksz cell geometry. The 
region of the cell in which the director field is close to homeotropic anchoring (9 = 0) increases 
with |C|. 

The behaviour of extensile aligning materials in a HAN geometry is reported in Figs. [14] and 
[15] for smaller and larger values of ( respectively. The spontaneous flow is asymmetric. Initially 
there are oppositely flowing slabs of liquid crystals, which distort the director field by creating 
homogenously aligned region separated by thin regions of homeotropic ordering. These profiles 
are then supplanted by an asymmetric quasi-Poiseuille flow, which resembles the response of a 
purely passive HAN cell to a pressure difference driven flow [41]. At larger values of ( the director 
profile throughout is close to the one obtained for a Freedericksz cell, with only a highly distorted 
boundary layer to satisfy the homeotropic anchoring at the top plane (z — L). 
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C. Spontaneous flow in two dimensions 



Thus far, all simulations reported here were performed in a quasi- ID geometry, where transla- 
tional invariance is assumed along x and y. The same simplification is often employed in numerical 
studies of passive liquid crystals (see many examples in Ref. I32L as well as e.g. Refs. [25, 41] 
for rheological studies); moreoever, as shown above they allowed us to check detailed analytical 
predictions (calculated at director- field or EL level) in exactly this geometry [Q]. It is clearly im- 
portant and interesting to consider whether there are additional spontaneous flow instabilities in a 
higher dimensionality. With periodic boundary conditions such instabilities must spontaneously 
break the translational invariance in x and y; we limit our attention to this case, but note that 
confining cell walls might also play an important role. 

We next present 2D simulations (L z = 100, L y = 100, L x = 0) in which we again have two 
parallel plates, normal to z; translational invariance along x is maintained but periodic boundary 
conditions are used to allow breakdown of this along the flow direction, y. We initialised the 
simulation with the director field along the y direction except for points along the mid-plane z = 
L/2, in which there was an alternating tilt of ±10° in stripes (the width of the initial stripes did 
not affect the steady state reached at the end of the simulations). 

Fig. [16] shows results for a moderate value of the activity parameter £ (0.001), for which the 
liquid crystal enters the spontaneously flowing active phase. Spontaneous flow appears as a pair 
of convection rolls which lead to a splay-bend in-plane deformation of the director field profile. 
The order parameter is to a good approximation constant (q ~ 0.5) throughout the sample. The 
threshold at which the spontaneous flow appears is smaller than the one found in the quasi- ID 
simulation (for which with the same parameters ( c ~ 0.002, see above). This is due to the fact 
that along y effectively homeotropic anchoring conditions are seen, and the active phase is entered 
for a smaller value of ( in this geometry. Note that, since at onset of the convection rolls there are 
exactly two of these in the periodic cell, the details of the transition may now depend sensitively 
on the aspect ratio of the cell. 

As we go deeper into the active phase, the number of convection rolls is, at early times in 
the simulations, larger (Fig. [T7](al,bl)). These convection rolls then split up, and the flow field 
acquires an out-of-plane component (i.e. there is flow along the x direction). After this happens, a 
number of vortices form which lead to a complicated flow which is accompanied by the formation 
of defects (of topological strength ±1/2) in the director field profile. The simulation, followed in 
Fig. [LZl does not lead to a steady state. It would seem plausible that the corresponding trajectories 
in phase space may be chaotic, but we have not attempted to test this directly. Moreover, once a 
nonzero x velocity has been acquired, there is a strong possibility of breakdown of translational 
invariance in x; to explore this would require fully 3D simulations. Note however that in this 
regime the structural length scale of the flow appears small on the scale of the simulation cell and 
therefore might cease to be sensitive to its shape. 
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V. DISCUSSION AND CONCLUSIONS 



We have presented a hybrid lattice Boltzmann algorithm to solve the equations of motion of 
an active nematic liquid crystal. In our equations the orientational degrees of freedom are char- 
acterised by a tensorial order parameter. This renders our algorithm general enough to deal - in 
principle - with non-homogeneous, flow-induced or paranematic ordering, as well as with topo- 
logical defects. The model we analyse is equivalent to the one proposed in Ref. 111]. 

Our main results are the following. First, we have explicitly mapped our model onto the one 
considered in Ref. |Q] in the limiting case of a uniaxial liquid crystal with a spatially uniform and 
time independent magnitude of ordering. This is useful when comparing the different approaches 
which are now being proposed to study the physics of active materials. 

Second, we found a spontaneously flowing phase (active phase) for a wide range of values 
for the activity parameter ( in a quasi- ID geometry where the director field is constrained to lie 
along a common direction along both confining plates. (A second activity parameter, A, merely 
renormalizes the equilibrium parameters of the passive material.) Our simulations confirm the 
location of the phase transition from passive to active phase found via a linear stability analysis 
in Ref. [7], but show that, for a wide range of parameters within the active phase, even very 
close to the boundary, the spontaneous flow profile has a quite different symmetry from the one 
predicted by that analysis. Instead of a sinusoidal flow with a node at the midplane, flow-aligning 
and flow-tumbling liquid crystals display a quasi-Poiseuille flow and a "boundary layer" type flow 
respectively. (Both flow profiles are bistable.) 

Our numerical method can readily probe, for the first time, the hydrodynamic behaviour of 
active materials deep in the active phase, where we gave evidence of a spontaneously banded 
flow for the flow-aligning case. Far from the phase boundary, there are multiple (initial condition 
dependent) solutions, and the system displays hysteresis. 

Third, if conflicting (HAN-type) anchoring conditions are applied at the confining plates, spon- 
taneous flow occurs for any values of the activity parameter (, however small. Finally, we per- 
formed two-dimensional simulations, with periodic boundary conditions along the y direction and 
planar anchoring along that direction on both confining plates. These suggest that there are addi- 
tional instabilities in a quasi-2D geometry. Moreover, at high activity levels, there can also be a 
spontanous flow also in the x direction in this geometry. 

These results demonstrate a remarkable richness in the steady-state hydrodynamic behaviour of 
active nematic materials, even in the absence of exernal drive such as an imposed shear flow. (As 
such, they have no counterpart in the physics of passive nematics.) Our hybrid lattice Boltzmann 
methodology, which combines LB for momentum with finite difference methods for the order 
parameter tensor Q a p, offers a robust and efficient method for probing these effects. It can equally 
well handle transient phenomena, some of which we explored above, and can readily be modified 
to allow for imposed flow. 

Our algorithm can be generalized in several ways. For instance, an additional order parameter 
equation, describing the time evolution of a polar vector field, can be considered with little more 
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effort. This would allow a full 2D study of polar active nematics [60 with a variable degree of 
ordering. Similarly, chiral active liquid crystals can be straightforwardly treated [42D, for instance 
to model concentrated actomyosin solutions. Actin fibers in very concentrated solutions undergo 
a nematic to cholesteric transition; another candidate for an active chiral liquid crystal might be 
a solution of DNA fragments interacting with polymerases or other motors [12]. Also, it would 
be of interest to use the present algorithm to characterise the rheological properties and map out 
the flow curves of an active liquid crystal under imposed shear. We shall report on such work in 
future publications. We also hope to report soon on fully three-dimensional simulations of active 
materials, along the lines pioneered for passive nematics in Il40ll . 

We acknowledge EPSRC for support, and are grateful to L. Tubiana for a critical reading of 
this work. 
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APPENDIX A: COMPARISON OF HYBRID WITH CONVENTIONAL LB CODES 



In Fig. [T8] we show the director and velocity dynamics at z = L/A (in the geometry of Fig. 



lybrid algorithm discussed in this paper 
40(1 for passive nematics and in [20] for 



[3]) and in the mid-plane respectively, computed via the 
and via a full LB algorithm (as described in Refs. ||39l . 
the active case). The agreement proves the validity of our hybrid approach. Note that two full 
LB algorithms are benchmarked against the hybrid code. In one case the double gradient term 
is entered as a constraint in the second moment, in the other its derivative is entered as a body 
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force (this second procedure guarantees that no spurious velocities are found in steady state, see 
e.g. Ref. |43D). It can be seen that the LB treatment with the double gradient terms entered in the 
second moment constraint leads to a small deviation at intermediate times. This we interpret as a 
discretisation error, as this method in 2D is known (for conventional i.e. passive liquid crystals) to 
lead to discretization errors causing small spurious velocities even in the steady state [14311 . 



APPENDIX B: "ERICKSEN-LESLIE" LIMIT OF THE ORDER PARAMETER EVOLUTION 
EQUATION 

In this Appendix we map the order parameter evolution equation used in this work, Eq. [51 
onto the analogous equation used in Ref. [7], by taking the limit of a uniaxial liquid crystal with 
spatially uniform and temporally constant magnitude of ordering q. In this way we will recover 
Eqs.[28]and[29l 

To this end let us first write the Q evolution equation © for H. This gives, formally, 

rH= {d t + u- V)Q-S(W,Q) -AQ. (Bl) 
By considering the uniaxial expression for Q (see Eq. ITTT) we obtain 

THfo = (dtq^pn^ ^d t q + (w 7 <9 7 g) npn^ ^ (w 7 <9 7 g) 5^ 

+q (d t np) n M + qn p {d t n,j) + q (u^np) n M + qnp (u^n^) 
6 2 

-Xqripn^ + \q^- + -£(q - l)D Plt 

—£q (Dj^njUfj, + npn^D Jfl ) - q (VL^n^n^ - npnyQ^) 

+2q£n f m ll Tr(QW) - ^(q - l)Tr(QW). (B2) 

As can be easily checked, one can substitute W with D in (|B2I) . As we have assumed that q does 
not depend on t and f, we obtain: 

THfa = qin^Np + npNp) - q^D^n^n^ + npn 7 D^) - Xq {npn^ - 

2 1 2 

+ g(<? - l)t D Pit + ^inpn^D^n^ + -q(l - q^Sp^D^n^ (B3) 

where Np, are co-rotational derivatives defined as, 

Np = dtnp + u^d^np + f2^ 7 n 7 

= d t np + u^np - {u x n)^ (B4) 

and u = V x u/2. In order to write the evolution equation (IB3I) in a form that resembles the one 
introduced in [7] we note first that, by the chain rule, 

bT bT dQ a p 



h 



Sriu 5Q a p dn 



H a pq(np5 afl + n a 5p^) 

qinpEp^ + n a H a ^) = 2q{npHp lx ). (B5) 
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If we now multiply (on the left) both members of Eq. (IB3I) by rip and we use the constraint 
n /3 n l3 = 1 we obtain after some algebra, 

Thj2q = q N, i --(q + 2)£n 7 D^ - ^Xqn^ (B6) 

where we have omitted terms 0(n 3 ). Clearly, if A = 0, Eq. (IB6I) reduces to the usual Ericksen- 
Leslie equation for the director field, namely Ii32ll 

hf, = jiN^ + j2n a D a ^ (B7) 

where 

2q 2 

71 = -jr, (B8) 

72 = -§(g + 2)e (B9) 
If, on the other hand, the active term A^Owe have 

2q 2 M+ 3 q 

Note that terms proportional to drop out of the equations in this mapping. Indeed they con- 
tribute a component of the molecular field parallel to n M , which would tend to increase the mag- 
nitude of the director field q. This is prevented by the Lagrange multiplier which appears in the 
vectorial "Ericksen-Leslie" model (to maintain constant q). As a result such terms simply change 
the relationship between the Lagrange multiplier and the magnitude of order and not the structure 
of the director field equation. By comparing Eq. (IB 101) with Eq. (3) of fl7|] (there Dn^/ Dt = N^) 
we then obtain the relations listed in Eqs. (l28l) and ( T2"9l ) in the text. 



iV M = — K + (BIO) 



APPENDIX C: "ERICKSEN-LESLIE" LIMIT OF THE NAVIER-STOKES EQUATION 

In this Appendix we work out the details of the mapping between the Navier-Stokes equation in 
our tensorial model in the uniaxial limit of constant q, and the momentum balance equation used 
in the "Ericksen-Leslie" version of Ref. [7], which was reported in Section III B in the text. To 
this end, we need to write the total stress tensor H af3 = ]j^ stve _|_ Y[^ lve in terms of the molecular 
and director fields, and respectively, which are used in director field based models. As in 
Appendix B we write Q a p in uniaxial form i.e Q = q(P — 1/3) where P a /3 = n a np. Note that 
P 2 = P and Tr(P) = 1 and recall that the Ericksen-Leslie expression for the total stress is: 

<rf^ = ain a npn^n p D^ p + a 4 D a( s + a^npn^D^ (CI) 
+a.Qn a n^D^ + a 2 n (3 N a + a 3 n a Np. 

We first consider the anti-symmetric part of the passive stress tensor in the tensorial model, namely: 

T a/3 = Q H - H Q 

= g(P-H-H-P). (C2) 
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Multiplying to the left the expression (IB3I) for H ai by P ay = n a n y and to the right by n 7 np, gives, 
after some algebra 



^ T af3 = gr(n a n 7 if 7j a — H ai n^np) 

q 2 {n a Np - N a n p ) - —{q + 2)(n a n 7 D j/3 - D^n^rip) 

Eq. (IC3I) may now be compared to the antisymmetric part of Eq. (IC2I) . to give 



(C3) 



2q 2 

«3 - «2 = -p- = 7i> (C4) 

"6 - "5 = = 72 (C5) 

where the equalities with 71, 72 come from comparison with Eq. ( IB9I ). We may slightly rewrite 
the antisymmetric term (IC3I) in a form that is closer to the one used in [Q]. This can be done by 
substituting the expression for iV M written in terms of the molecular field 

K_ ll n<jD ( C6) 
7i 7i 

into (IC3I) . This gives 

rr a/3 = — (n a hp - h a n p ) + q 2 — {n a D aa n p - n a n a D a(3 ) - — (q + 2)(n a n a D a(3 - D aa n a np). 
7i 7i <3 

(C7) 

Hence the expression for r Q/ g simplifies to 

q 2 

r a/3 = — - {n a h p - h a n p ) (C8) 
1 71 

which is the antisymmetric term in the director field treatment of Ref. [0] (see eq. (2) of [0]). 

We now turn to the symmetric part of the total stress tensor (excluding the active contribution 
and the double gradient term): 

o"a/3 = - PoS af3 + 2£(<3 Q/ 3 + -5 a p)Q lt H le (C9) 
The active contribution is: 



n^ ive = -C^^ + C§<W. (cio) 



Note that the double gradient term term —d a Q~ jV S9 q is analogous to the director field term 
~da n v S Q~ F n > which is not included in Eq. (IC2I) hence not considered hereafter. 
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By using Eq. (|B3I) for H, after some algebra, one obtains the complete expression for a a p as 

-^(q + 2)(n f3 N a + n a Np) + ^ 



&af3 = ~^p(9 + 2 )( n f3 N a + n a Np) + — (4 - q) (D^n^np + n a n 1 D 1 i 3 ) 



4 2 

gp(5 - iffAtf + 7^<? 2 £ 2 (V - 4g - 3)n a npD lv n u n„ 



2 



+ ^(4-7g-8g 2 + 8g 3 )5 a/3 D 7 ^n 



'7 



The first term of the right hand side of Eq. (IC11I) can be usefully rewritten (for comparison 
with the equation in QTU) by using (IC6I) to write A 7 ^ in terms of h^. 



<?£/3(<? + 2){n p N a + n^iVg) = -q£/3(q + 2) f n/j— - n p —n a D aa + - n a — ) 

V 7i 7i 7i 7i / 

(n^/ia + n a ^) — (n/3n a D aa + n a n a D al3 ){C\\) 



2 

where in the last line we have used Eq. (f29b . 

The Navier Stokes equation in the Stokes regime is 

ridp(d a U{] + d p u a ) = 2dpr]D al3 = -dp(H. a p). (C12) 

—H a p can equivalently be rewritten as 

- Il Q/3 = --^p (nph a + n a hp) + -^r(q ~ l f {p m n a np + n a n a D af3 ) 

4 2 

- gp(9 - l) 2 £ 2 A*/3 - 7^g 2 £ 2 (V - 4g - 3)n a npD^ v n v n^ 

- —(4 - 7g - 8g 2 + Sq^^D^n^ 

+ (qn a np - (-5 a p 
q 2 

- — (n a hp - h a np) . (C13) 
171 

If A = ( = 0, i.e. for passive liquid crystals, Eq. (IC11I) gives the symmetric part of the Beris- 
Edwards stress (ignoring the distortion stress) and this, together with Eq. (|C3I) gives the Leslie 
coefficients which are listed in Section III B (Eqs. l3lT[3~6l) . 

In Eq. (IC13I ) the term proportional to D a/3 may be added to the left hand side in Eq. IC 121 to 
renormalise the apparent viscosity, while the rest of it may be rewritten as 



— {nph a + njip) + — (<? - l) 2 (D aa n a np + n a n a D afi ) 

2 

- 7^P<? 2 £ 2 ( 4 <? 2 - 4g - 3)n ol npD^ u n v n 1 

- —(4:-7q-8q 2 + 8q 3 )5 al3 D ll ,n u n 7 
+ Cqn a nf3 - c|5 Q/3 

- ^ {n a h/3 - haUp) (C14) 
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where for the last term we have used relation (IB9I ). By comparing our equation with the one in 
we then get Eqs. (1301) . (1291) in the text. 



23 



(a) (b) 




-1.6 - 



7 Y 

FIG. 1: Plot of the 7 dependence of vel\ panels in (a) and (b) have \ = 0.7, 0.5 respectively. Within each 
panel different curves refer to different activity levels A (see legend). Note that for £ = 0.5, flow tumbling 
< 1) is exnected throughout the nematic nhase. Points A and R renresent numerical examples 
described below. 
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FIG. 2: Phase boundary for L = 49 in the (A, Q plane for 7 = 3.0, Tt = 1 and £ = 0.7. Four points found 
numerically from our HLB simulations are also shown (filled circles). 
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FIG. 3: Geometry used for the calculations described in the text. The active gel is sandwiched between two 
infinite plates, parallel to the xy plane, lying at z = and z = L. We consider (a) normal anchoring and (b) 
conflicting anchoring. (The latter would correspond to a hybrid aligned nematic (HAN) cell for a passive 
liquid crystal material.) 
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FIG. 4: Time evolution of the components of the polarization field n y (upper) and n z (lower), at z = L/4. 
Parameters are L = 100, ( = 0.005, A = , 7 = 3, t/ = 2.5 and £ = 0.7 (flow aligning regime). At the 
bounding plates, the field is strongly anchored along the y direction (homogeneous anchoring). 
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FIG. 5: Profiles of director orientation angle (a) and velocity field (b; in lattice units) at steady state for 
different values of £ in a flow-aligning active liquid crystal sample with L = 100 (other parameters as spec- 
ified in the text). Solid, dashed and dot-dashed curves correspond to £ = 0.003, 0.005, 0.01 respectively. 
The transition to the active phase occurs at ( = ( c ~ 0.002. The flow is bistable: reversing the sign of 9 
and u y together creates an alternative steady-state solution. 




FIG. 6: Profiles of director orientation (a) and velocity field (b, lattice units) at steady state for different 
values of ( in a flow-aligning active liquid crystal sample with L = 100 (other parameters as specified in 
the text). Solid, dashed and dot-dashed curves correspond to £ = 0.02, 0.04, 0.08 respectively. All solutions 
are bistable (see text). 
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FIG. 7: Profiles of director orientation (a) and velocity field (b, lattice units) at steady state for different 
values of ( in a flow-aligning active liquid crystal sample with L = 400 (other parameters as specified in the 
text). Solid, dashed, dot-dashed and dotted lines correspond to ( = 0.001, 0.002, 0.003, 0.01 respectively. 
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FIG. 8: Time evolution of the components of the polarization field at z = L/2 for the flow tumbling case 
(£ = 0.5). Other parameters are L = 100, ( = -0.0025, A = 0, A = 0.1, and K = 0.04; the transition as 
predicted by Eq. @2]is at £ = £* ~ —0.0022. The director field is strongly anchored along the y direction 
(homogeneous anchoring). 
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FIG. 9: Polarization angle (a) and velocity field (b) profiles for flow-tumbling active liquid crystals, with 
C=-0.003 (solid black line), -0.004 (dashed red line), -0.005 (dot-dashed green line), and -0.006 (dotted blue 
line). The transition between the passive and the active phase is attained at £ = Q ~ —0.002 (see also Eq. 
S2]>. 
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FIG. 10: Polarisation angle (a) and velocity (b) profiles for flow-tumbling active liquid crystals deep in the 
active phase. Curves correspond to £=-0.008 (solid black line), -0.01 (dashed red line), -0.02 (dot-dashed 
green line), and -0.03 (dotted blue line). 
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(a) (b) 




FIG. 1 1 : Profiles of director orientation (a) and velocity (b) for two different steady state solutions found for 
contractile tumbling liquid crystals in the active phase (£ = —0.003) in the geometry of Fig. [3^), starting 
with two different initial conditions. In (a) the solid and the dot-dashed line refer to the two different 
polarisation angles, while the long dashed line refers to the <fi angle between the projection of the director 
angle onto the xy plane and the positive x axis. Initial conditions are given in the text. 
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FIG. 12: Profiles of director orientation (a) and velocity (b) for two different steady state solutions found 
for contractile tumbling liquid crystals in the active phase (£ = —0.006) in the geometry of Fig. [3^). Initial 
conditions are given in the text. 
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FIG. 13: Profiles of director orientation (a) and velocity (b) for flow-aligning contractile active liquid crys- 
tals in a HAN geometry. Curves correspond to £ = —0.001 (solid black line), —0.0005 (dashed red line), 
-0.003 (dot-dashed blue line). 




FIG. 14: Profiles of director orientation (a) and velocity (b) for flow-aligning extensile active liquid crystals 
in a HAN geometry. Curves correspond to £ = —0.001 (solid black line), —0.0005 (dashed red line), 
-0.003 (dot-dashed blue line). 
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FIG. 15: Profiles for director orientation (a) and velocity (b) for flow-aligning extensile active liquid crystals 
in a HAN geometry. Curves correspond to ( = —0.001 (solid black line), —0.0005 (dashed red line), 
-0.003 (dot-dashed blue line). 




FIG. 16: Maps of velocity field (a) and director field (b) in steady state for an active aligning liquid crystal 
with ( = 0.001 (extensile), simulated on a two-dimensional L = 100 x L = 100 grid. 
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FIG. 17: Maps of velocity field (al-a3) and director field (bl-b3) for an active aligning liquid crystal with 
( = 0.01 (extensile), simulated on a two-dimensional L = 100 x L = 100 grid. The three rows correspond 
to the configurations after 10 4 , 3 x 10 4 , 10 5 lattice Boltzmann steps respectively. 
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FIG. 18: Time evolution of Q yy at z = L/4 (a) and of u y in the mid-plane (b), as predicted by our hybrid 
LB treatment (solid black lines), and by two types of full LB treatment (dashed red lines, with the double 
gradient terms entered in the first moment constraint, to avoid spurious velocities at equilibrium; and dot- 
dashed blue lines, with the double gradient term entered in the second moment constraint). 
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